Density-matrix approach for an interacting polariton system 
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Using Lindblad approach we develop a general formalism for theoretical description of a spatially 
inhomogenous bosonic system with dissipation provided by the interaction of bosons with a phonon 
bath. We apply our results to model the dynamics of an interacting one dimensional polariton 
system in real space and time, analyzing in detail the role of polariton-polariton and polariton- 
phonon interactions. 



I. INTRODUCTION 

A semiconductor microcavity is a photonic structure 
designed to enhance the light-matter interaction. In a 
planar microcavity photons are confined between two 
mirrors and resonantly interact with the excitonic tran- 
sition of a two-dimensional semiconductor quantum well 
(QW). If the quality factor of the cavity is sufficiently 
high, it is possible to achieve the regime of strong cou- 
pling between the cavity photon and QW exciton. In 
this case the elementary excitations in the system, which 
are called cavity polaritons, have a hybrid, half-light, 
half-matter nature. The peculiar properties of polaritons 
make them a unique laboratory for studying of various 
collective phenomena interesting from the point of view of 
basic physics, which range from polariton BEG super- 
fluidity y and Josephson effect [J] to polariton-mediated 
superconductivity [5|. 

Besides fundamental interest, quantum microcavities 
in the strong coupling regime can be used for a vari- 
ety of optoelectronic applications @. Recently, it was 
proposed that the peculiar spin structure and possibil- 
ity to achieve lateral confinement of polaritons |7| opens 
a way for creation of optical analogs of spintronic com- 
ponents (so-called spinoptronic devices), based on trans- 
port of cavity polaritons in real space. In this context, the 
analysis of onc-dimensional (ID) polariton transport is of 
particular importance Q, as ID polariton channels are 
fundamental building blocks of such future spinoptronic 
devices as p olariton neurons Q and polariton integrated 
circuits [10[. 

On the theoretical side, transport properties of 
exciton-polaritons in real space have not yet been studied 



in detail. Early works jlll - ll3j treated the case of polari- 



tons interacting with external potentials only neglecting 
both polariton-polariton and polariton-phonon interac- 
tions. However, in the most interesting regime of polari- 
ton BEG neither of the two interactions can in principle 
be neglected. Indeed, coupling of the polaritons with a 
reservoir of acoustic phonons leads to thermalization of 
the polariton subsystem, which is dramatically speeded 
up by polariton- polariton interactions which are known 
to be responsible for overcoming the so-called "bottle- 



neck effect" [ij]. Besides, in the regime of polariton 



BEG polariton-polariton interactions are responsible for 
the onset of superfluidity [2[. 

Gurrently, there are two ways of describing a sys- 
tem of interacting polaritons. Assuming full coherence 
of the polariton system, polariton-polariton interactions 



can be accounted for using nonlinear Gross- Pitaevskii 
(GP) equation, satisfactorily describing the dynamics of 
inhomogeneous polariton droplets in real space and time 
[l7| . The approach, however, does not include in- 
teractions with a phonon bath, responsible for thermal- 
ization of the system and leading to dephasing. In the 
opposite limit, when the polariton system is supposed to 
be fully incoherent, its dynamics can be described us- 
ing a system of semiclassical Boltzmann equations 
2l| . which provides information about time dependence 
of the occupation numbers in reciprocal space but fails 
to describe real space dynamics in the inhomogeneous 
system. 

Recently, there appeared theoretical attempts to com- 
bine the two mentioned approaches introducing dissipa- 
tion terms into GP equation in a phenomenological way 
[2^ . To our mind, however, these attempts, although 
being interesting and leading to a rich phenomenology 
23I . 24 1, are not fully satisfactory, as they lack any mi- 
croscopic justification. In the present paper we give an 
alternative way of describing the dynamics of an inhomo- 
geneous polariton system in real space and time account- 
ing for dissipation effects. Our consideration is based on 
Lindblad approach for density matrix dynamics. We use 
our results for modeling of the propagation of polariton 
droplets in ID channels. It should be noted, that the 
method we develop is rather general and can in prin- 
ciple be applied to any system of interacting bosons in 
contact with a phonon reservoir (for example, indirect 
excitons [25|). We have previously applied this formal- 
ism to model Josephson oscillations between two spatially 
separated condensates, but in that context the spatial de- 
pendence was trivial [26|. 



II. FORMALISM 

We describe the state of the system (polaritons plus 
phonons) by its density matrix p, for which we apply 
Born approximation factorizing it into the phonon part 
which is supposed to be time-independent and corre- 
sponds to the thermal distribution of acoustic phonons 

Pph = exp |—/3 i/ph I and the polariton part ppoi whose 

time dependence should be determined, p = pph ® Ppoi- 
Our aim is to find dynamic equations for the time evo- 
lution of the single-particle polariton density matrix in 
real space. 

p(r,r',t) =rr(^t(r,^)^(r',t)p| ^ {i^\v,t)i^{v' ,t)) (1) 
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where '0^(r, t), ^/'(r, t) are operators of the polariton field, 
and the trace is performed by all degrees of freedom of the 
system. The particularly interesting quantities are the 
diagonal matrix elements which give the density of the 
polariton field in real space and time n{r,t) = p{r,r,t). 
In our consideration we neglect the spin of the cavity po- 
laritons for simplicity, as our main goal here is to account 
for dissipative dynamics coming from interaction with 
phonons, and spin degree of freedom is not expected to 
introduce any qualitatively new physics from this point of 
view. However, that introduction of spin into the model 
is straightforward and corresponding work is currently 
underway. It should also be noted that our model is to 
some extent similar to one proposed in Ref. 1271 . Differ- 
ently from that paper, however, we assume that all de- 
coherence in the system comes from the interaction with 
acoustic phonons which is neglected in Ref. |27| and thus 
do not perform the separation of the polariton ensem- 
ble into coherent low energy part and incoherent high- 
energy reservoir. 

It is convenient to go to reciprocal space, making a 
Fourier transform of the one-particle density matrix. 



p(k,k',t) = {27rf/L'^ J e''-^''-^''''^p{r,r',t)drdr' = (2) 

= Tr{a+ak+q} = {a+a^') 

where d is the dimensionality of the system (d = 2 for 
non-confined polaritons, d — 1 for the polariton channel), 
L is its linear size, a^, Ok are creation and annihilation 
operators of the polaritons with momentum k. Note, that 
we have chosen the prefactor in a Fourier transform in 
such a way, that the values of p(k, k', t) are dimensionless, 
and diagonal matrix elements give occupation numbers 
of the states in discretized reciprocal space. Knowing 
the density matrix in reciprocal space, we can find the 
density matrix in real space straightforwardly applying 
inverse Fourier transform. 

The total Hamiltonian of the system can be repre- 
sented as a sum of two parts. 



of the polariton- polariton interaction constant coming 
from Hopfield coefficients 15|. We do this approxima- 



H — Hi + H2 



(3) 



where the first term Hi describes the "coherent" part 
of the evolution, corresponding to free polariton propa- 
gation and polariton-polariton interactions and the sec- 
ond term Hphon corresponds to the dissipative interac- 
tion with acoustic phonons. The two terms affect the 
polariton density matrix in a qualitatively different way. 



tion because the goal of the manuscript is to present a 
novel formalism for description of the relaxation effects 
and not detailed modeling of a particular experiment, 
and we want to keep our formalism as simple as possible. 
Besides, this approximation is widely used in current de- 
scription of polaritonic systems based on modifications 
of the Gross-Pitaevskii equations 22l-l24l|. However, for 
modeling of realistic experiments the p-dependence of 
the interaction constant can easily be introduced into the 
equations. 

The effect of Hi on the evolution of the density matrix 
is described by the Liouvillc-von Neumann equation. 



i {dtp) 



(1) 



[Hi;p] 



(5) 



which after use of mean field approximation leads to 
the following dynamic equations for the elements of the 
single-particle density matrix in reciprocal space (see Ap- 
pendix I for details of the derivation): 

{dtp{K k)}^'^ = 2C/ ^ Im {p(ki, ki - p)p(k, k + p)} 

(6) 

{atp(k,k')}^'^ = i{E^ - i?(.)p(k,k')+ (7) 
+ ^C/ 5^ p(ki, ki - p) [p(k - p, k') - p(k, k' + p)] 

ki,p 

These expressions represent an analog of the Gross- 
Pitaevskii equation written for the density matrix. 



B. Scattering with acoustic phonons 

Polariton-phonon scattering corresponds to the inter- 
action of the quantum polariton system with a classical 
phonon reservoir. It is of dissipative nature, and thus 
straightforward application of the Liouville-von Neu- 
mann equation is impossible. Introduction of dissipation 
into quantum systems is an old problem, for which there 
is no single well established solution. In the domain of 
quantum optics, however, there are standard methods 
based on the Master Equation techniques [l^- In the 
following we give a brief outline of this approach applied 
to a dissipative polariton system. 

The Hamiltonian of interaction of polaritons with 
acoustic phonons in Dirac representation can be repre- 
sented as 



A. Polariton-polariton interactions 

The part of the evolution corresponding to Hi is given 
by the following expression 

i^i = ^ -Bkajok + ^ ^ a+ a^^aki+pak2-p (4) 

k ki,k2,p 

where i?k is the energy dispersion of the polaritons, U 
is the matrix element of the polariton-polariton interac- 
tions. In the current paper we neglect the p-dependence 



H2{t) = H-{t)+H+{t) 



(8) 



i(Bk+q--Ek)i 



'k+q 



ak(6qe-'"-* + 6iqe* 



k,q 



where Ok are operators for polaritons, 6q operators for 
phonons, and Wq are dispersion relations for po- 
laritons and acoustic phonons respectively, -D(q) is the 
polariton-phonon coupling constant. In the last equality 
we separated the terms H^ where a phonon is created, 
containing the operators 6"*", from the terms H^ in which 
it is destroyed, containing operators b. 
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Now, one can consider a hypothetical situation when 
polariton-polariton interactions are absent, and ah re- 
distribution of the polaritons in reciprocal space is due 
to the scattering with a thermal reservoir of acoustic 
phonons. One can rewrite the Liouville-von Neumann 
equation in an integro-differential form and apply the so 
called Markovian approximation, corresponding to the 
situation of fast phase memory loss (see Ref. l28| for the 
details and discussion of limits of validity of the approx- 
imation) 



[dtp) 



(2) 



dt'[H2mH2{t');pm= (9) 

= 5^E [2 (H+pH- + H-p H+) - 

- {H+H- + H-H+) p- p {H+H- + H-H+)] 



where the coefficient 6^e denotes energy conservation 
and has dimensionality of inverse energy and in the cal- 
culation taken to be equal to the broadening of the po- 
lariton state [l^ . For time evolution of the mean value of 
any arbitrary operator (A) = Tr{pA) due to scattering 
with phonons one thus has (derivation of this formula is 
represented in Appendix II): 

{dt{A)y'^ = 6ae {{[H-; [A; H+]]) + {[H+: [A; ff-]])) . 

(10) 

Putting A = a^ai^i in this equation we get the contri- 
butions to the dynamic equations for the elements of 
the single-particle density matrix coming from polariton- 
phonon interaction (see Appendix III): 



{a*p(k,k)}(^) = 



(11) 



J2 2W{q') [p(k+q',k+q')(</ + l)(p(k,k) + 1) 

q',Sk<-Bk+q' 



-p(k,k)n^^(p(k+q',k+q') + i; 
J2 2I^(q') [P(k+q', k+q>^';, (p(k, k) + 1) 

q',Bk>-Ek+<,, 



p(k,k)«'' , + l)(p(k+q',k+q') + 1) 



and 



(12) 



ph 



X< J2 W^(q') [P(k+q',k+q') 
+ ^(q') [-P(k+q', k+q') - < - 1 

q',£;k>-Ek+q' 

+ W)[p(k' + q',k' + q')-< 

q',-E;,<-Bk'+q' 



+ J2 ^(q') [-pO^' + q'' + q') - < - 1" 

where the transition rates are given by Vr(q) — 
-D^(q)i5A_Eq ^-iid nj^'' denote the occupation numbers of 
phonons with wavevector q given by Bose distribution. 



Eq. pT|) is nothing but a standard Boltzmann equa- 
tion for the phonon-assisted polariton relaxation, while 
Eq. (fT2|) describes the decay of the off-diagonal matrix 
elements of a single-particle density matrix due to inter- 
action with classical phonon reservoir. Together these 
equations thus describe thermalization of the polariton 
system. 

To account for the effects of free polariton propaga- 
tion, polariton-polariton and polariton-phonon interac- 
tions one should combine together expressions ([5]), ([7|), 
(|TT]) and (|T2|) . After finding the single-particle density 
matrix in reciprocal space by solving the corresponding 
dynamic equations, one can determine the dynamic of the 
system in real space simply performing a Fourier trans- 
formation by k and k' variables. 



III. RESULTS AND DISCUSSION 

Our formalism is suitable for the description of both 2D 
polaritons and polaritons confined within ID channels. 
In this paper we present the results of numerical mod- 
eling for the latter case only, as solving of the dynamic 
equations for 2D case requires the use of supercomputing 
facilities. 

The results of modeling are shown on Figures 1-3. We 
consider a 2 pui wide polariton channel in GaAs mi- 
crocavity with Rabi splitting 15 meV at temperature 
T = 1 K. The polaritons are created by a short coherent 
localized laser pulse. We account for the finite lifetime of 
cavity polaritons r = 5 ps adding the term — p(k, k')/r 
into the dynamic equations. 

Figure 1 shows the dynamics of the polariton system 
in reciprocal space and demonstrates the roles of the 
polariton-phonon and polariton-polariton interactions. If 
only the former are included, the system demonstrates 
the bottleneck effect shown on Fig. la. Due to the energy 
relaxation coming from the interactions with phonons, 
polaritons have a tendency to move towards the ground 
state in fc-space. However, this process is dramatically 
slowed down in the inflection region of the polariton dis- 
persion, where polariton-phonon interaction becomes in- 
efficient. Consequently, there is no remarkable increase of 
the population of fc = state [l3l ■ The bottleneck effect 
can be overcome by the polariton-polariton interactions, 
as shown on Fig. lb. One sees that in this case the par- 
ticles accumulate quickly in fc = state. At the same 
time, the second maximum of the polariton distibution 
appears at higher k due to the energy conservative na- 
ture of the polariton-polariton interactions (analogically 
to the formation of the idler mode in polariton paramet- 
ric oscillator [30j). 

Our results for the dynamics of the polariton distri- 
bution in reciprocal space are in good agreement with 
those obtained by using Boltzmann equations. In addi- 
tion, our approach allows consideration of the dynamics 
of the dissipative polariton system in real space. This is 
illustrated on Figures 2 and 3. 

Figure 2 shows the effect of the various types of interac- 
tions on the real space dynamics of the localized polariton 
wavepackage. We compare the cases of the ballistic prop- 
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FIG. 3: Polariton distribution in real space at times t = 0.5 
ps (inset) and t — 1.5 ps after creating the package by a lo- 
calized laser pulse centered around fc = 3 • 10^ m~^. Scat- 
tering on phonons and polaritons are both accounted for. 
The temperatures are: 1 K (red/solid), 4 K (green/dash- 
dot), 8 K (blue/dashed), 15 K (magenta/dotted) and 20 K 
(black/solid). 



Time tps) 



Wavevector (m l ! 



FIG. 1: Time evolution of polariton distribution in fc-space. 
The vertical plane denotes k — state, a) Only polariton- 
phonon scattering is accounted for. The relaxation towards 
fc = state is blocked due to the bottleneck effect, b) Both 
polariton-phonon and polariton-polariton scattering are ac- 
counted for. The maximum of the polariton concentration is 
developed at fc = 0, signifying the overcoming of the bottle- 
neck effect. Another maximum appears at higher fc, due to 
the energy-conserving nature of the polariton-polariton inter- 
actions. 
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FIG. 2: Polariton distribution in real space at times t — 0.5 
ps (inset) and t = 1.5 ps after creating the package by a local- 
ized laser pulse centered around fc = 0. Black solid line corre- 
sponds to ballistic propagation, red dashed line - to polariton- 
phonon interactions, blue dotted line - to polariton-polariton 
interactions and green dashed line - to both polariton-phonon 
and polariton-polariton interactions. 



agation with those where only polariton-phonon interac- 
tions are included, only polariton-polariton interactions 
are included and both polariton-polariton and polariton- 
phonon interactions are included. As one sees, the dy- 
namics are very different for these four cases. Polariton- 
polariton interactions lead to splitting of the wavepack- 
age into two soliton- like peaks, which is in good qualita- 
tive agreement with the results given by Gross-Pitaevskii 
equation [l7| . On the other hand, polariton-phonon in- 
teractions lead to damping of the package, contributing 
to the recovering of the homogenious distribution of the 
polaritons in real space as it is expected from the classical 
diffusion equation. 

Naturally, the effect of the phonon damping strongly 
depends on temperature, as shown on Fig. 3. One sees, 
that at low temperatures the propagating package is split 
into two due to the polariton-polariton interactions. In- 
creasing the temperature smoothes down the polariton 
distibution and at T = 20K one has just a single peak. 



IV. CONCLUSION 

In conclusion, we developed a formalism for the de- 
scription of the dissipative dynamics of an inhomoge- 
neous polariton system in real space and time account- 
ing for polariton-polariton interactions and polariton- 
phonon scattering. The formalism was applied for 
numerical modeling of the propagation of a polariton 
droplet in a ID channel. Our approach can also be used 
for modeling the dissipative dynamics of other bosonic 
systems (e.g. indirect excitons). 
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APPENDIX I: DERIVATION OF KINETIC 
EQUATIONS FOR POLARITON- POLARITON 
SCATTERING 

The evaluation of the dynamic equations of (aj^flk) and 
(a]jj.ak+q) using Eq. ([Sj and applying the mean field ap- 

I 



Tr {a+Ok+q [p; H.p^p]} = U ^ p{k2, ka - p)Tr {a+Ok+q [p; a+ ak^+p] } = (13) 

ki ,k2 ,p 

= [/ ^ p{k2, k2 - p) i-pik, k+q+p) + p(k-p, k+q)] 
k2,P 



TT{a^ai,[p;Hp-p]} ^ /?(k2, ka - p)Tr {a+ak [p; aj^aki+p] } = (14) 

ki ,k2 ,p 

= X! P(k2,k2-p)Tr{p[a|[^aki+paJak-ai^akaJ^aki+p]} = -2iL/ ^ Im {p(k2, ka - p)p(k, k+p)} 

ki,k2,p k2,p 

From this one straightforwardly gets Eqs. El [T] 



APPENDIX II: DERIVATION OF EXPRESSION delta-function omitted) if dynamics of a density matrix 
FOR DYNAMICS OF MEAN VALUES IN BORN- is given by Eq. U 
MARKOV APPROXIMATION 

Now, consider the evolution of a mean value of any 
arbitrarty operator A, (A) = Tr{pA) (energy conserving 

I 



dt{A) = Tr {2H+pH-A + 2H- pH+ A - {H+R- + R- H+) pA - p {H+R- + R- R+) A) = (15) 
= Tr [p {2R-AR+ + 2R+AR- - A {R+R- + R-R+) - {R+R- + R-R+) A)] 

where we used the property of the invariance of the trace as regards cyclic permutations of operators. The latter 
expression can be simplified as 

2R-AR+ + 2R+AR- - A {R+R- + R-R+) - {R+R- + R-R+) A = (16) 
= 2R-{[A; R+] + R+A) + 2R+{[A- R-] + R- A) - A {R+R- + R-R+) - {R+R- + R-R+) A = 
= 2R-[A:H+\ +2R+[A;R-] - R+[A;R-] - [A;R+]R- -R-[A;R+] - [A:+R-]R+ = 
^[R-;[A;R+]] + [R+;[A;R-]] 

I 

where we used the following property of commutators: 
[A;BC] = B[A;C] + [A;B]C. Thus 

lR-;[A;R+]]+ = -[[A;R+] + ;R+]= (18) 
= [[R- ; A+]; R+] = [R+ :[A; R-]] 

dt{A) = Tr {p[R-; [A; R+]]) +Tr {p[R+; [A; R-]]) (17) 



For the important case of the Hermitian operator A 
corresponding to a physical observable one has 



at(A) = 2Re{Tr {p[R- ■,[A; R+]])} 



(19) 
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This formula can be applied for calculation of occupa- 
tion numbers. 

APPENDIX III: DERIVATION OF DYNAMIC 
EQUATIONS WITH ACOUSTIC PHONONS 

To get explicit expressions for dynamics of p(k, k, t) let 
us consider a simple case when only states k, k+q arc 
present. We have two cases: 

a) ^^(k + q) > E(k). In this case, leaving energy- 
conserving terms only one gets 

H+ = Da+a^+^b+ (20) 
H- = Da+ ai^bq (21) 



The application of EqUU] gives the following results: 



S-^dtp{K k) = S-^dtUi, = Tr {p[H-; [a+at; H+]]) + 

(22) 

Tr {p[H+; [a+a^; R-]]) = 2Tr {p[H+; [a+a^; R-]]) 



After some straightforward algebra one gets 



[a+ak;H ] = D[a^ai^; a^^^a]^bq} = -Da+_^q&q[a+ak; = -Da+^^aj^b^, (23) 

[a+ak;-ff"]] = --D^[a+ak+q&q ; a+_^qakfoq] = 

= (a^+qak+q(akak + l)ib^K + 1) - («^+q«k+q + l)c4akb+bqj 

and 

dtp{k, k) = 2SD^ [p{k + q, k + q)« + l)(p(k, k) + 1) - p{k, k)nP'^(p(k + q, k + q) + 1)] (24) 

I 



This is nothing but ordinary Boltzmann equation, ac- 
counting for transition from state k+q to state k ac- 
companied by the emission of the phonon (spontaneous 
or stimulated) and transition from state k to state k+q 



due to the phonon absorption. 

b) i?(k + q) < £^(k). Treating this case in a similar way 
we find: 



^(^(k, k) = 2SD^ p{k + q, k 



q)n^_'l(p(k,k) 



l)-(p(k + q,k + q) 



l)(n^_"q + l)p(k,k) 



(25) 



r 



Again we obtained Boltzmann equation, but differently 
to the previous case the transition from state k+q to 
state k goes with absorption of the phonon and from 
state k to state k+q with phonon emission. Performing 
summation over all reciprocal space one gets Eq. ([TT|). 

Now let us consider the dynamics of the off-diagonal 
part of the density matrix p(k, k' ,t). 



S-^dtp{k,k') = 
--Tr{p[H+;[a+a^,;H-]])+Tr {p[H- 



(26) 



akflk';^^]] 



Here we should consider different orderings of the en- 



ergies corresponding to the states k, k' and state k + q. 
As an example, let us consider the case when < 
i?k+q, £'k' < ^-k+q- Leaving energy-conserving terms 
only one gets 



= Datau+nb, 



qOq +ak,ak+qOq, 



bt 



(27) 

= i'akflk+q&q + flk'aJ+qV (^8) 
where q' = k + q — k'. In this case the first term in 
Eq. dMl) gives: 



[afcOkS-ff ] 
[H+-[a+ai,r,H-]] 



= -D[aJak';aJ+qak&q + ak'aJ+qV] = ^-^V^k'aJ+q 
= D'^[ai,+cia^b+ + flk+qOk'^q'; -^q«k+q«k'] = 

= D'^ [a^a^'ia^+qa^^+q + l)^'q6q - {b+bq + l)a+_^qak+qa+ak| 



(29) 
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and 



= [au+q^k+qa+ak>{b+bq, + 1) - (ak+qOk+q = l)akak'&q'V} 



(30) 



Finally for the case iJk < ^'k+q, E\^' < -Ek+q one obtains equation in the form 



dtp{k, k') = <5Z?2 (p(k + q, k + q)« + 1) - (p(k + q, k + q) + l^') p(k, k') + 
p(k + q, k + q)(<'/ + 1) - (p(k + q, k + q) + p(k, k') 



(31) 
(32) 



The same procedure is applied for all other cases. Per- 
forming again summation over all reciprocal space one 
gets Eq. [m 
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